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I. PROTECT OVERVIEW 


I.A. Scientific Objectives In order to understand the physical and chemical processes which 
produce the tenuous planetary and planetary satellite (upper) atmospheres through interactions 
with their particle, field, and radiation environs, it is necessary analyze remotely observed and 
spacecraft data with physically meaningful models. With this in mind, we have undertaken a 
coupled program of theoretical modeling and complementary data analysis regarding the global 
distributions of neutral and ionized gases in, and escape from, tenuous planetary atmospheres of 
Io and Europa. The theoretical models developed will have further applications to other tenuous 
atmospheres such as that of Pluto and Titan or the upper atmospheres of the terrestrial planets. 

I.B. Significance and Relevance This project addresses prime areas of study central to current 
and future NASA missions and for the goals of NASA’s overall program in the study and 
exploration of the solar system and its origin. In the published report of COMPLEX, the 
Committee on Planetary and Lunar Exploration, of the Space Studies Board, National Research 
Council, comets and the Jupiter system were identified in the highest priority group. 

Io and Europa, their atmospheres, and their interactions with Jupiter’s plasma torus, are 
prominent in current and future NASA science goals and missions. A unique combination of 
modeling tools which have been, and continue to be, developed at the University of Michigan 
have already been applied to understand Galileo particle and fields measurements obtained 
during its encounters with Io and Europa. Both the Galileo mission and the Galileo Europa 
mission and remote measurements with the Hubble Space Telescope will continue to be sources 
of new results, as data already taken is reduced and released, and as upcoming encounters with Io 
and Europa yield new important insight into these atmospheres. NASA is developing plans for 
the Europa Orbiter mission for studying the icy surface and mantle below the icy surface, which 
may hold the clues to the thin extended atmosphere. 

I.C. Strategy and Approach Our overall strategy is and has been to develop and apply 
sophisticated models of the type normally used in purely theoretical studies for the purpose of 
directly comparing models with measurements in order to gain a deeper understanding of the 
physical state and composition of the atmosphere in question. Typical first-principles theoretical 



investigations produce many detailed calculations to explore the type and range of physical 
phenomena that might or should occur, but are often only loosely constrained by observed 
characteristics and are often limited to physically overly simplistic assumptions. On the other 
hand, typical measurements, whether Earth-based and near-Earth based remote observations, or 
in situ spacecraft measurements, are often limited to fairly simple model analyses (scale-height 
distributions, uniform columns, constant temperatures, simple imposed uniform magnetic fields, 
etc.). The Pi’s approach throughout his studies with collaborators dealing both with comets and 
satellite atmospheres has been to use physically meaningful models and compare directly with 
data. For this reason we have developed the core modeling tools, and either analyze published 
measurements, or establish various collaborations with observers. The subjects of study covered 
here are unified by this overarching strategy. 

II. Progress Report 

Work under this grant has lead to publication of two papers on the MHD modeling of the 
interaction of the plasma torus with Io (Combi et al. 1998; Kabin et al. 2000), two papers on the 
MHD modeling of the interaction of the plasma torus with Europa (Kabin et al. 1999; Liu et al. 
2000), a review paper on a large set of MHD model results, (Kabin et al. 2000), and the writing 
and submission of a refereed review paper for a special AGU monograph on the comparative 
aspects of tenuous cometary and planetary satellite atmospheres/ionospheres and their 
interactions with their outer particles and fields environs— either planetary magnetosphere or 
solar wind (Combi, Gombosi and Kabin 2002). In addition, some of the ion particle-pushing 
algorithms included in our kinetic DSMC model were described in a paper on the transport of 
particles from the impact site of comet Shoemaker-Levy 9 on Jupiter to the conjugate 
hemisphere (Bauske, Combi and Clarke 2000). Dr. Konstantin Kabin, obtained his Ph.D. in 2000 
(Kabin 2000) with much of his support from this grant. 

Refereed publications describing work supported by this grant are: 

1) Combi, M.R., K. Kabin, T.I. Gombosi, D.L. DeZeeuw, K.G. Powell. lo’s Plasma 
Environment during the Galileo Flyby: Global Three-Dimensional MHD Modeling with 
Adaptive Mesh Refinement. J. Geophys. Res. 103 (A5), 9071-9081 1998. 



2) Bauske, R., M.R. Combi, J.T. Clarke. Analysis of Mid-Latitude Auroral Emissions 
Observed During the Impact of Comet Shoemaker-Levy 9 With Jupiter . Icarus 142 , 106- 
115, 1999. 

3) Kabin, K., M.R. Combi, T.l. Gombosi, A.F. Nagy, D.L. DeZeeuw, K.G. Powell. On 
Europa’s Magneospheric Interaction: An MHD Simulation of the E4 Flyby. J. Geophys. 
Res. 104 , 19983-19992, 1999. 

4) Kabin, K., K.C. Hansen, T.L Gombosi, M.R. Combi, T.J. Linde, D.L. DeZeeuw, C.P.T. 
Groth, K.G. Powell, A.F. Nagy. Global MHD Simulations of Space Plasma 
Environments: Heliosphere, Comets, Magnetospheres of Planets and Satellites. 
Astrophys. Space Sci. 274 , 407-421, 2000. 

5) Liu, Y., A.F. Nagy, K. Kabin, M.R. Combi, D.L DeZeeuw, T.L Gombosi, K.G. Powell. 
Two-Species 3D MHD Simulation of Europa's Interaction with Jupiter's Magnetosphere. 
Geophys. Res. Lett. 27, 1791-1794, 2000. 

6) Kabin, K„ M.R. Combi, T.L Gombosi, K.C. Hansen, K.G. Powell. Io's Magnetospheric 
Interaction: An MHD Model with Day -Night Asymmetry. Planet. Space Sci. 49, 337-343, 

2001. 

7) Combi, M.R., T.l. Gombosi, K. Kabin. Plasma Flow Past Cometary and Planetary 
Satellite Atmospheres, in Comparative Planetary Aeronomy, AGU Monograph Series ( in 
press), 2002. 

Preliminary calculations using the kinetic DSMC model, given here, have been described 
at a number of conferences (Bauske and Combi 1998a&b, 1999). Although we would not 
describe our recent publication output from the current grant as inadequate, we do acknowledge 
that publication of a full paper documenting the model and giving first results has been hampered 
by a number of unfortunate non-technical problems. These include loss of a Co-I to industry 
because of uncertainty in continued funding and a delayed start of a graduate student because of 
personal circumstances. A new student, Mr. Valeriy Tenishev, has come to Michigan with a M.S. 
degree, having significant experience with DSMC methods and related programming. He has 
been coming-up-to-speed in learning the DSMC code. In the remainder of the current grant year 
we will complete a study running all neutral species in the DSMC model and adopting the ion 
flow from our published MHD models, with the goal of submitting a paper for publication. 



Impact of This Work . The combination of modeling tools at our disposal provides us with a 
unique capability to address all the important components of a coupled system of neutral gas 
molecules, radicals and atoms, plasma ions and electrons and electromagnetic fields. Most 
modeling efforts concentrate on one aspect of the problem, and do that part extremely well. 
However, most so-called “self-consistent” treatments of one aspect of the problem, no matter 
how well-done, are forced to make rough assumptions about one of the other fundamentally 
important parts of the problem. 

Copies of all of the refereed publications are attached as an appendix to this report. A 
summary of the kinetic ion and neutral DSMC modeling effort (Bauske & Combi, 1998a&b, 
1999) follows. 

Kinetic Ion-Neutral DSMC Model . Important processes in thin tenuous atmospheres, such as that 
of those around Jupiter’s satellites, like Io, and the upper atmospheres and ionospheres of the 
terrestrial planets, work on various energy and rate scales. In order to merge the many local 
processes into a global three-dimensional picture, a model is needed, which can include these 
processes on a small spatial scale and which can diversify the different energy scales. At the 
same time there are disparate time scales for various collision and gyration rates. This is difficult 
with a continuum (hydrodynamic) model alone. In order to model acceleration processes an 
electrodynamic or PIC type model is need, but then there are difficulties to include the 
chemistry. Our hybrid approach is a way out of this situation. 

We have developed a 3D neutral and ion multi-species kinetic model, embedded in the 
multi-scale Cartesian grid of a 3D MHD code. The MHD part of this code has already been used 
in simulations of the mass-loaded flow of Jupiter's corotating magnetospheric plasma past Io 
(Combi et al., 1998) and other planetary applications (Bauske et al., 1998a, b; Kabin et al. 1999, 
2000 & 2001; Liu et al. 2000). The DSMC part is specifically targeted towards analysis of 
kinetic processes as expected in tenuous atmospheres, from the ground outward to a distance of 
few times the body's radii. 

The model couples the Direct Simulation Monte Carlo (DSMC) method (Bird, 1994) 
with a magnetohydrodynamics (MHD) (Powell et al., 1999). DSMC has a long history of 
application in areas as diverse as study of hydrodynamic shock structure, calculation of low- 
density, hypersonic flows around spacecraft and simulations of plasma reactors for 



microelectronics manufacturing. Originally, the method was used to simulate the transition 
regime, where the mean free path of particles is too large for continuum hydrodynamics to be 
applicable. Because collisions are important, free molecular simulations are not appropriate. 
Therefore, individual particles are simulated as they move around within a grid, colliding with 
other particles and with solid objects. Macroscopic properties, such as density and temperature 
are computed by appropriate averaging of particle masses, locations, velocities, and internal 
energies. Momentum and energy exchanges with the surface allow for chemical reactions and 
sputtering effects. DSMC is based on the “rarefied-gas” assumption that over a short time 
interval or 'step' the molecular motion and the intermolecular collisions are uncoupled and 
therefore can be calculated independently. Molecules are moved over the distances appropriate 
for this time step, followed by the calculation of a representative set of collisions. The time step 
is small compared to the mean collision time and, the results are independent of its actual value. 

Our model is divided into inner and outer domains. As shown in Figure 1, inner domain 
may contain a whole planetary obstacle (Io) and a fraction of its surrounding atmosphere / 
exosphere / plasma environment, e.g., a cube of a few Io radii. Both types of domains are 
simulated using MHD to obtain the magnetic field. In the inner domain, we embed the multi- 
species DSMC model for detailed simulations of neutral and plasma chemistry and energetics 
and surface interactions. An adaptively refined unstructured grid gives high resolution in all 
areas of interest. Each domain contains cells of various sizes to match specific local 
requirements. The MHD part can have high resolution at a planetary bow shock, while the 
DSMC part is able to simulate the outflow of hot gases near a volcanic source. At domain 
boundaries, the MHD flow entering a DSMC domain is converted to individual particles by 
sampling from the given velocity distribution. The particle-tracking algorithm accounts for 
gravity, centripetal, Coriolis and electromagnetic forces. Particles are tracked through both 
domains, collisions, however, are only calculated within the DSMC domain. 

Our approach is unique in several ways: (1) It is a new method to incorporate a 
microphysical perspective into the continuum calculations of magnetohydrodynamics. (2) It 
differs considerably from pure electromagnetic approaches like the Particle-In-Cell (PIC) 
method, because PIC codes usually neglect particle-particle and particle-boundary collisions, or 
include Coulomb collisions without chemistry (Winske and Omidi, 1996). (3) Many applications 
of the DSMC method today are in dilute, but collisionally dominated particle regimes where the 




Figure 1 . Hybrid model for lo. 


particles have nearly thermal collision energies and the ionization ratio is very low. This is not 
the case for the tenuous atmospheres. Our model can also explore the gap between the collision 
dominated and the collision free regime in partially ionized plasmas. We extend the range of the 
DSMC method towards the higher particle energies usually encountered in planetary ionospheric 
and magnetospheric plasma environments (hundreds of eV/amu). 

Lutisan (1995) compared several variants of the DSMC method and found the No Time 
Counter technique (NTC) of Bird (1994) and the Null-Collision technique (Koura, 1986 & 1989) 
to be among the fastest and most reliable ones. We chose the NTC method for its simplicity and 
its closeness to classic kinetic theory. In a DSMC cell of volume V c each simulated molecule 
represents F n real molecules. The probability P of collisions between two simulated molecules 
during the time interval At is equal to the ratio of the volume swept out by their total cross- 
sections o T moving at the relative speed g between them to the cell volume, 

P “ F n o T gAt / V c . 






The average number of simulated molecules is N = nVJF n where n is the density of the real 
molecules. P is generally a very small quantity. Therefore it is inefficient to calculate the full set 
of collisions by selecting all of the N(N-l)/2 simulation pairs and computing the collisions with 
this probability. The procedure can be made more efficient if only a fraction of the possible 
collision pairs is included and the resultant probability increased by dividing by this fraction. 
Efficiency is maximized if the fraction is such that the maximum probability becomes unity. 
Bird's NTC method for collisions of species p molecules with molecules of species q is that 

XN p N q F n (o r g)™ At/V e 

pairs are selected and the collisions are computed with probability 

p P_ o T g 

(o r g)Z 

where N q is the average number of species q simulation particles in the cell and o T g is the 
product of cross section and relative velocity of the selected pair. Note that the rate of collisions 
is not affected by the value of (o T g )™' . This parameter is initially set to a large but reasonable 

value for this species combination and stored for each cell with the provision for it to be 
automatically updated if a larger value is encountered. After certain numbers of collisions, it is 
reset in order to account for time dependent changes in the species population. 

In contrast to Bird, we allow a particle to undergo several changes of its internal energy 
per time step, thus allowing for fewer particles per cell and larger time steps where necessary. 
We added functionality beyond the basic algorithms suggested by Bird (1994) to account for 
chemical reactions and collisions with large inelastic cross sections. Apart from higher energies 
than in usual DSMC calculations, we also face the difficulty to simulate regions with orders of 
magnitude differences in density and also the inclusion of a model chemistry, which contains a 
variety of different species including trace species. Our basic idea for solving these problems 
numerically is to use weighted species groups and a regularly spaced subgrid within each cell of 
the adaptive grid structure. Usually every particle in a simulation cell represents a fixed number 
of real particles. This leads to problems if one species is underrepresented since a large number 
of samples is necessary to get reasonable statistics for this species. A weighting method (e.g., see 
Serikov, 1991; Miller and Combi, 1994) artificially changes the species value in order to allow 
for more simulation particles and to keep the overall number computationally feasible. 



We developed a weighting 
scheme consistent with the NTC 
method, which allows for trace species 
and conserves the overall momentum 
and energy balance on average. The 
number of selections for collisions is 
recalculated every time a species type- 
changing collision (reaction, 
ionization, etc.) occurs and the 
collision iteration proceeds as long as 
there are selections remaining. Our 
NTC variant employs a data structure, 
which allows fast addition and 
deletion of particles without needing 
to re-index all particles after 
collisions. It determines the number of selections and the acceptance of the collisional 
momentum, energy, and type-change for the higher weighted particle. Weighting is used to keep 
computationally reasonable numbers of simulation particles in different density regimes roughly 
equal by comparing particle weights to a cell weight (Combi, 1996). Particles with higher 
weights are replicated, and those with lower weights are deleted with a probability set by the 
ratio of the particle weight to the cell weight. 

We have tested and included methods (Boyd, 1993 & 1994; Wysong and Wadsworth, 
1998) to calculate the probability of internal energy redistribution during a collision based on the 
species combination, the individual collision energy, and the internal states of the participating 
particles. The particle-based approach is more reliable at high energies and low densities than 
methods that use constant or temperature-dependent probabilities. Reliable DSMC methods are 
required to reach equipartition of internal and translation energies under equilibrium conditions. 
This requirement is fulfilled in our methods, if constants of the underlying molecular models are 
properly calibrated. Figure 2 shows an example for rotational and vibrational relaxation of 0 2 . 

Our hybrid model specified for Io, includes species O, S, 0 2 , S 2 , SO, S0 2 , S0 3 , and the 
corresponding ions. We account for 0 + ( 4 S), 0 + ( 2 D), 0 + ( 2 P), S + ( 4 S), S + ( 2 D), S + ( 2 P) as separate 
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Figure 2. Rotational and vibrational relaxation of 0 2 . Starting 
with a translational temperature of 4700° and without internal 
energy a particle distribution maintained in a box with 
reflecting walls relaxes to a common temperature indicating a 
detailed balance in all internal states. Ttrans, Tvib, and Trot are 
the translational, vibrational and rotational temperatures, 
respectively. The DSMC model includes various methods to 
calculate these temperatures (Bird 1994, Borgnakke and Larsen 
1975, and Boyd 1994, Hansen 1991). 




species for photoionization and 
take care of the excess energies, 
which contribute to the electron 
energies for ionization, to the 
particle energies for dissociation, 
and to UV heating. Energy 
dependent solar photo rates and 
the solar photon fluxes are taken 
from Huebner et al. (1992). The 
plasma is currently assumed to be 
quasi-neutral, and as a first 
approximation the electron 
temperature is given by the scaled 
ion temperature— an improvement 
we plan to address. 

Figure 3 shows 
intermediate results of a 
simulation for Io, where the 
evolution of an S0 2 sublimation 
and volcanic atmosphere without 
influences from the Io torus is 
studied. We simulated a region of 
30 Io radii starting with a vacuum. 
Here we show S0 2 densities in a 
selected region around Io and on 
lo's surface. After 1000 iterations using local time steps for each cell, the outflow of gas from the 
dayside of the planet and at several volcanoes on the day- and night-side almost filled the 
simulation region, but a state of equilibrium is not yet established. However, from comparisons 
with previous steps we know that changes in the maximum density between succeeding 
iterations are small and decreasing close to the surface, while the remaining empty cells fill with 



Figure 3. SO z densities in a sublimation and volcanic atmosphere 


around Io. The sun is to the left and a number of volcanic sources is 
distributed on lo's surface. Here, we show intermediate results of a 
DSMC, which started with a vacuum. After 1000 iterations the 
average densities in the lower altitudes of the atmosphere are 
established (compared to preceding results). The results for density 
and velocity are similar to the hydrodyanamic results of Moreno et 
al., however the temperatures at high altitude are higher. This is 
similar to the results of Austin and Goldstein (2000) who found 
significant kinetic effects not too far above the surface. The linear 
cell size of our adaptive grid is 0.05 Io radii and each cell contains a 
subgrid with 64 cells to resolve the lower atmosphere. 





simulation particles, and average values of variables such as cell densities temperature and flow 
velocities are established. 

Finally, the PI would like to gratefully acknowledge the Planetary Atmospheres program 
and the current and past program managers for their support for this work. 
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